%5-dec use esk-hua instead
%4-sep-2007. This code reads Julia and ACe data sets
%pick days with cont. common measurements
% 
load c:\manoj\projects\plasma\HUA_PIU1.mat
load c:\manoj\geomag\indices\aplist.mat;
load c:\manoj\projects\ace\ESKLER fday_esk_ler_ts_filt FILT_LER FILT_ESK;
esk_fday = fday_esk_ler_ts_filt;
%FILT_ESK = (LER_TS-nanmean(LER_TS)+ESK_TS-nanmean(ESK_TS))./2;
FILT_ESK = (FILT_LER+FILT_ESK)/2;
%clear fday_esk_ler_ts LER_TS;
clear fday_esk_ler_ts_filt FILT_LER;
ESK_floor_fday = floor(esk_fday);
load c:\manoj\geomag\indices\aplist.mat;
mean_delta = (60/(24*3600));


Hd = chebi11_filter();
N_seg = 1;
N_data = 0;
mjd_date = datenum(2000,1,1);
for i = 1: size(HUAPIU,1),
 %   if sum(isnan(HUAPIU(i,:))) <= 60, 
    L = ace_all(:,1) >= HUA_PIU_fday(i)+0.40 & ace_all(:,1) <= HUA_PIU_fday(i)+0.8;
    %L = ace_all(:,1) >= HUA_PIU_fday(i)+0.555 & ace_all(:,1) <= HUA_PIU_fday(i)+0.8333;
    if sum(L) > 0,
        ace_time = ace_all(L,1)+(20)/(60*24);%Try to advance 20 ace time minutes
        %ace_time = ace_all(L,1);
        ace_data = ace_all(L,2);
        mean_imf_bz = nanmean(ace_gsm_bz(L));
        bz_data = ace_gsm_bz(L);
%
        le_sl = length(ace_data);
        L = isnan(ace_data);
        if sum(L) < 10 & sum(L) ~= 0,
            ace_data = interp1(ace_time(~L),ace_data(~L),ace_time);
        end;
            fprintf('ACE fday %d ACE missing %d ', HUA_PIU_fday(i), sum(L));

        L = isnan(HUAPIU(i,:));
        if sum(L) < 50 & sum(L) ~= 0,
            HUAPIU(i,:) = interp1(magnetic_frac(~L),HUAPIU(i,~L),magnetic_frac);
        end;
        
%
    fprintf(' HUA fday %d HUA missing %d\n', HUA_PIU_fday(i), sum(L));
   % if mean(diff(ace_time)) <= 0.0029,
   if abs(mean(diff(ace_time))-0.0035) <= 2.7778e-005,
        HUAH = interp1(HUA_PIU_fday(i)+magnetic_frac, HUAPIU(i,:), ace_time);
        
% UGLY

        HUAH(1) = HUAH(2);
% UGLY OVER        
        
        ACED = ace_data;
         HUAH1 = filter(Hd,HUAH);
         ACED1 = filter(Hd,ACED);
       
%         subplot(311);
        L = fday_ap >= HUA_PIU_fday(i) & fday_ap < HUA_PIU_fday(i)+1;
        mean_ap = mean(ap(L));
        if sum(isnan(ACED)) == 0 & sum(isnan(HUAH)) == 0 & mean_ap >= 20,% & mean_imf_bz < 0 ,
            HUAH_SEG(N_seg:N_seg+le_sl-1) = HUAH;
            ACE_SEG(N_seg:N_seg+le_sl-1) = ACED;
            TIME_SEG(N_seg:N_seg+le_sl-1) = ace_time;
            N_seg = N_seg+le_sl;
            N_data = N_data+1;
            data_length(N_data) = length(HUAH);
%                
          end;
   else,
%        fprintf('Day %d has some missing time stamp\n', Julia_W(i).fday);
    end;
    end;
%end;
end;

% %[Cxx,F] = mscohere(JULI_SEG,ACE_SEG,hanning(72),0,72,0.0033333);
[Cxx,F] = mscohere(HUAH_SEG,ACE_SEG,hanning(le_sl),0,le_sl,0.0033333);
%  
figure1 = figure;

axes('Parent',figure1,'XTick',[10 20 30 50 70 100 200 300 500 1000],...
    'XScale','log',...
    'XMinorTick','on');
set(gca,'FontSize',16);

box('on');
hold('all');
xlabel('period in minutes');
ylabel('coherence');
hold on;
 semilogx((1./(60*F)),Cxx,'r','LineWidth',2);
 axis([10 1000 0 1]);

% 
% 
alpha = 0.95;
ci = 1- (1-alpha)^(1/(N_data-1));
aa=axis;
h=line([aa(1) aa(2)],[ci,ci],'LineStyle','-.','color','r');
